
% HHVR_prediction_2t.m is called by HHVR_2016_Master_code.m
 
 function  [percent_correct_ML  percent_correct_EM  ] = HHVR_prediction_2t(  parameters_ML,  parameters_EM, z )  ;

% First, take the probabilities using CLogit parameters
   
  
global DATA choices   y_short X     individuals
 
s=size(DATA,1)/choices ;
num1=zeros(s,1);num2=zeros(s,1);num3=zeros(s,1);num4=zeros(s,1);num5=zeros(s,1);
prob_1=zeros(s,1);prob_2=zeros(s,1);prob_3=zeros(s,1);prob_4=zeros(s,1);prob_5=zeros(s,1);
den= zeros(s,1); den_x1= zeros(s,1);den_x2= zeros(s,1);den_x3=zeros(s,1);
score = zeros(s,size(parameters_ML,2));
 
 for i=1:(size(DATA,1))/choices
  
     num1 (i) =  exp(X(i*choices-4,1)*parameters_ML(5) + X(i*choices-4,2)*parameters_ML(6)     ) ;
     num2 (i) =  exp(X(i*choices-3,1)*parameters_ML(5) + X(i*choices-3,2)*parameters_ML(6)   + 1*parameters_ML(1)  ) ;
     num3 (i) =  exp(X(i*choices-2,1)*parameters_ML(5) + X(i*choices-2,2)*parameters_ML(6)   + 1*parameters_ML(2)  ) ;
     num4 (i) =  exp(X(i*choices-1,1)*parameters_ML(5) + X(i*choices-1,2)*parameters_ML(6)   + 1*parameters_ML(3)  ) ;
     num5 (i) =  exp(X(i*choices  ,1)*parameters_ML(5) + X(i*choices  ,2)*parameters_ML(6)   + 1*parameters_ML(4)  ) ;
     
       den(i)   =  num1 (i)  +  num2 (i) + num3 (i) + num4 (i)  +  num5 (i)   ;
       prob_1(i)     =  num1 (i) /  den(i)  ;
       prob_2(i)     =  num2 (i) /  den(i)  ;
       prob_3(i)     =  num3 (i) /  den(i)  ;
       prob_4(i)     =  num4 (i) /  den(i)  ;
       prob_5(i)     =  num5 (i) /  den(i)  ;

 end
 
 prob_ML = [prob_1 prob_2 prob_3 prob_4 prob_5];
  prediction_ML = zeros(size(prob_ML,1), size(prob_ML,2));
 for i = 1 : size(prob_ML,1)
     max_ = max([prob_ML(i,1) prob_ML(i,2) prob_ML(i,3) prob_ML(i,4) prob_ML(i,5)]);
     for j = 1 : choices
         if prob_ML(i,j) >= max_
             prediction_ML(i,j)=1;
         end
     end
 end
 clear max_
 
%checking
sum_ = sum(prediction_ML) ; sum_ = sum(sum_) ; 
if sum_ ~= individuals
    display 'something is wrong predicting CLogit'
    pause
end

prediction_short_ML = zeros(individuals,1);
 for i = 1 : size(prob_ML,1)
     for j = 1 : choices
         if prediction_ML(i,j) >= 1
             prediction_short_ML(i,1)=j;
         end
     end
 end

comparing_ML = [y_short prediction_short_ML];

num_correct_ML_ =   zeros(individuals,1);
 for i = 1 : size(prob_ML,1)
      if comparing_ML(i,1)==comparing_ML(i,2)
        num_correct_ML_(i,1)=1;
      end
 end
percent_correct_ML_=sum(num_correct_ML_)  % /individuals; 
size(percent_correct_ML_,1)
percent_correct_ML = percent_correct_ML_ /size(num_correct_ML_,1) ;
 
 %          Assign to each observation a type with probability 1, according to z .
 %          That is, take the 3 probabilities that we have, and assume that
 %          this observation fails in the category for which the
 %          probability is largest
 
 
 PI = mean(z);
 type_= zeros(s, 2);
 for i = 1 : size(prob_ML,1)
     max_ = max([z(i,1) z(i,2)    ]);
     for j = 1 : 2
         if z(i,j) >= max_
             type_(i,j)=1;
         end
     end
 end
  
 clear max_
%checking
sum_ = sum(type_) ; sum_ = sum(sum_) ; 
if sum_ ~= s
    display 'something is wrong predicting heterogeneous types'
    pause
end
type = zeros(s,1);
 for i = 1 : size(type,1)
     for j = 1 : 2
         if type_(i,j) >= 1
             type(i,1)=j;
         end
     end
 end 
  
 
 s=size(DATA,1)/choices ;
num1=zeros(s,1);num2=zeros(s,1);num3=zeros(s,1);num4=zeros(s,1);num5=zeros(s,1);
prob_1=zeros(s,1);prob_2=zeros(s,1);prob_3=zeros(s,1);prob_4=zeros(s,1);prob_5=zeros(s,1);
den= zeros(s,1); den_x1= zeros(s,1);den_x2= zeros(s,1);den_x3=zeros(s,1);
score = zeros(s,size(parameters_ML,2));
size(score);
 
for i=1:(size(DATA,1))/choices
  
     if type(i)==1
     num1 (i) =  exp(X(i*choices-4,1)*parameters_EM(5)     ) ;
     num2 (i) =  exp(X(i*choices-3,1)*parameters_EM(5)   + 1*parameters_EM(1)  ) ;
     num3 (i) =  exp(X(i*choices-2,1)*parameters_EM(5)   + 1*parameters_EM(2)  ) ;
     num4 (i) =  exp(X(i*choices-1,1)*parameters_EM(5)   + 1*parameters_EM(3)  ) ;
     num5 (i) =  exp(X(i*choices  ,1)*parameters_EM(5)   + 1*parameters_EM(4)  ) ;
     
       den(i)   =  num1 (i)  +  num2 (i) + num3 (i) + num4 (i)  +  num5 (i)   ;
       prob_1(i)     =  num1 (i) /  den(i)  ;
       prob_2(i)     =  num2 (i) /  den(i)  ;
       prob_3(i)     =  num3 (i) /  den(i)  ;
       prob_4(i)     =  num4 (i) /  den(i)  ;
       prob_5(i)     =  num5 (i) /  den(i)  ;
     elseif type(i)==2
     num1 (i) =  exp(X(i*choices-4,2)*parameters_EM(6)     ) ;
     num2 (i) =  exp(X(i*choices-3,2)*parameters_EM(6)   + 1*parameters_EM(1)  ) ;
     num3 (i) =  exp(X(i*choices-2,2)*parameters_EM(6)   + 1*parameters_EM(2)  ) ;
     num4 (i) =  exp(X(i*choices-1,2)*parameters_EM(6)   + 1*parameters_EM(3)  ) ;
     num5 (i) =  exp(X(i*choices  ,2)*parameters_EM(6)   + 1*parameters_EM(4)  ) ;
     
       den(i)   =  num1 (i)  +  num2 (i) + num3 (i) + num4 (i)  +  num5 (i)   ;
       prob_1(i)     =  num1 (i) /  den(i)  ;
       prob_2(i)     =  num2 (i) /  den(i)  ;
       prob_3(i)     =  num3 (i) /  den(i)  ;
       prob_4(i)     =  num4 (i) /  den(i)  ;
       prob_5(i)     =  num5 (i) /  den(i)  ;
 
     else
         display 'not good'
     end
 end
prob_EM = [prob_1 prob_2 prob_3 prob_4 prob_5];
 
 
 
 prediction_EM = zeros(size(prob_EM,1), size(prob_EM,2));
 for i = 1 : size(prob_EM,1)
     max_ = max([prob_EM(i,1) prob_EM(i,2) prob_EM(i,3) prob_EM(i,4) prob_EM(i,5)]);
     for j = 1 : choices
         if prob_EM(i,j) >= max_
             prediction_EM(i,j)=1;
         end
     end
 end
 clear max_
 
%checking
sum_ = sum(prediction_EM) ;
sum_ = sum(sum_) 
if sum_ ~= individuals
    display 'something is wrong predicting CLogit'
    pause
end
 
prediction_short_EM = zeros(individuals,1);
 for i = 1 : size(prob_EM,1)
     for j = 1 : choices
         if prediction_EM(i,j) >= 1
             prediction_short_EM(i,1)=j;
         end
     end
 end
 
comparing_EM = [y_short prediction_short_EM];
 
num_correct_EM_ =   zeros(individuals,1);
 for i = 1 : size(prob_EM,1)
      if comparing_EM(i,1)==comparing_EM(i,2)
        num_correct_EM_(i,1)=1;
      end
 end
percent_correct_EM_=sum(num_correct_EM_);  % /individuals; 
size(percent_correct_EM_,1)
percent_correct_EM = percent_correct_EM_ /size(num_correct_EM_,1) 

